Stability and Dynamics of Crystals and Glasses of Motorized Particles 
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Many of the large structures of the ceU, such as the cytoskeleton, are assembled and maintained 
far from equilibrium. We study the stabilities of various structures for a simple model of such a 
far-from-equilibrium organized assembly in which spherical particles move under the influence of 
attached motors. From the variational solutions of the manybody master equation for Brownian 
motion with motorized kicking we obtain a closed equation for the order parameter of localization. 
Thus we obtain the transition criterion for localization and stability limits for the crystalline phase 
and frozen amorphous structures of motorized particles. The theory also allows an estimate of 
nonequilibrium effective temperatures characterizing the response and fluctuations of motorized 
crystals and glasses. 
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Assemblies of molecular sized particles are seldom far 
from equilibrium owing to the relative strength of the 
thermal buffeting inherent at this scale. As we con- 
sider assemblies of larger and larger particles the ther- 
mal forces become less capable of moving and reorganiz- 
ing such assemblies. At the size scale of biological cells, 
objects are not rearranged just by equilibrium thermal 
forces but are moved about by motors or by polymeriza- 
tion processes that use and dissipate chemical energy . 
What are the rules that govern the formation of peri- 
odically ordered or permanently organized assemblies at 
this scale? Does the far-from-equilibrium character of 
the fluctuating forces due to motors and polymer as- 
sembly change the relative stability of different colloidal 
phases? These problems are not unique for intracellular 
dynamics, but belong to an emerging family of nonequi- 
librium assemblyproblems ranging from driven parti- 
cles 0, swarms U, and jamming to microscopic 



pattern formation and mesoscopic self-organization 

Motivated by these considerations which may be rele- 
vant for the dynamics of the cytoskeleton P, Q as well 
as other far-from-equilibrium aggregation systems, we 
study a simple motorized version of the standard hard 
sphere fluid often used to model colloids. Both motors 
and nonequilibrium polymer assemblies can convert the 
chemical energy of high energy phosphate hydrolysis to 
mechanical motions which one would ordinarily think 
would "stir" and hence destabilize ordered structures. 
We will show these systems in some circumstances may 
have an enlarged range of stability relative to those with 
purely thermal motions. 

We adopt a stochastic description of the motions of 
a collection of motorized particles. The overdampcd 
Langevin dynamics is r,; = PDfi + ■q{t) + . Here r.i 
is the position of the ith particle, fi = —ViU is the 
mechanical force that comes from the usual potential 
U{{r}) := L/(ri,f2, . . . ,r„) = Y.{ij)^{^ij) among par- 
ticles. The random variable fj vanishes on average and 
is Gaussian with (vf{t)T^^^{t')\ = 2D6a0S,j6{t - t'). The 



motor term = J2q iq5{t—tq) is a time series of shot- 

noise-like kicks. Its properties depend on the underlying 
biochemical mechanism of the motors. The stochastic 
nature of the motors also leads to a master equation de- 
scription 1^, 1^ for the dynamics of the probability distri- 
bution function ^' of the particle configurations, 
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Here Lpp := DJ2i^i ■ (V» - is the Fokker- 

Planck operator. An integral operator LNE^{{r}) = 



J n,dr\ [A-^pj^^-jVl,({^'}) _ X^^,^^^p^vl/({f})j sum- 

marizes the nonequilibrium kicking effect of the motors. 

The motors are firmly built in the particles. They work 
by consuming chemical energy sources, like ATP. In a 
single chemical reaction event, the motor makes a power 
stroke (which induces a discrete conformation change) 
that moves the particle by a distance of £ in the direction 
n. Motor kicking can be modeled as a two-step stochas- 
tic process: step I, the energy source binds to the motor, 
and the step II, the reaction ensues and the resulting 
conformational change makes a power stroke. The rate 
of the first step ki depends on the energy source concen- 
tration while the rate of the second step /c2 depends on 
the coupling between the structural rearrangement and 
the external forces, = Kexp(s/3[C/(r) — U{r + £)]), i.e., 
motors slow down when they work against mechanical 
obstacles. Such slowing has been demonstrated in the 
case of microtubules s, the coupling strength, 

measures the relative location of the transition state for 
the power stroke step and ranges from to 1. At the 
limit s ^ 1 we have a susceptible motor and while s ^ 
corresponds to an adamant motor. We use these names 
in the sense that an adamant motor is not sensitive to its 
thermal- mechanical environment, so each power stroke 
uses and wastes a lot of energy; in contrast a susceptible 
motor saves energy running faster downhill (free energy) 
and slower uphill. 

We assumed that step II is the bottleneck of kinet- 
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ics, i.e., the overall rate fc w ^2 <C fci. To make our 
model suitable for a variety of situations, we specify 
different variables s and s' for the degree of suscepti- 
bility for going uphill and downhill respectively. Thus 
k = K exp{~ f3Q[U {f + £) — U{r)]) with a switch function 
G{x) := 'd{x)sx + Q{~x)s'x. Here Q is the Heaviside 
function. Thus the overall temporal statistics of the kicks 
is position-dependent Poisson distribution. 

The kicking direction n{t) fluctuates on the time scale 
r of the particle tumbling. We will study below explicitly 
two extremes: the isotropic kicking case when r is very 
small and the persistent kicking case when r is very large 
compared with k~^, i.e., each motor always kicks in a 
predefined direction. The direction of persistence will be 
assumed to vary randomly from particle to particle. 

To solve the dynamics of probability distributions, re- 
searchers often pose the problem as the solution of a 
variational problem. Due to the L^e part of eqn. (1), 
we can not perform the usual transformations of the left 
and right state vector to make L hermitian For this 
type of problems, nevertheless we can obtain the solution 
of the many-particle master equation using nonhermitian 
variational methods as described by Eyink 0, or us- 
ing the squared (therefore hermitianized) operator L^L 
(e.g., pl59of i). 

Eyink's nonhermitian variational formulation is sim- 
ilar to the Rayleigh-Ritz method in ordinary quantum 
mechanics but uses independent left and right state vec- 
tors. For isotropic kicking, we start with a Jastrow-likc 
trial function 

= exp |- - R^r] - PUi{r})^ (2) 

Similar to the quantum hard spheres [T^ . such a trial 
function avoids any singularities of L pp arising from the 
hardsphcre potentials Uij{rij) between particles i and j 
used in this study. For simplicity, we set ql := fi — Ri and 
a uniform := ^. The nonhermitian variational method 
implies for steady states that the second moment of q 
satisfies the moment closure [l5| equation 
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To effectively evaluate eqn. (3), we need to simplify the 
many-body integration J IliCpqi invol ving exp — /3[/({r}). 
Here we use cluster expansion 0, IVq to render the 
manybody Boltzmann factor a product of effective single 
body terms by averaging over the neighbors' fluctuations. 
Wc thus have e'^^ = H^e"''"' w U.e''^^^ Here the orig- 
inal Ui = '^j^i depends on the manybody configu- 
ration, while Ui depends on qi (and constant {i?}) only. 
In fact, we keep it to the harmonic order for consistency, 
i.e., Ui = P~^aqf. Here a is the effective spring constant 
from the mechanical feedback from neighbors, a depends 
on its neighbors' overall fluctuations controlled by a and 
their mean position {R}. In turn, the positions of the 



neighbors are controlled by the lattice spacing for crystals 
or radial distribution functions for glasses and ultimately 
by the nature of the structure and the particle density n. 
I.e., for fee lattice, we have acr = ctcr{ct, n) as the eigen- 
values of the Hessian matrix constructed from the effec- 
tive potential i J2j£,i n '^)- Here v{R; a) = ln{l+ 

ierf[(i? ~ l)V5]-ierf[(i? + l)Vfi] + i2Z^[e-"(«-i)'- 
g-Q(fi+i) jj. ^j^^ ^Yie sum of Rj is over the 12 nearest 
neighbor positions of the origin of a fee lattice (with 
lattice spacing (^)^). For glasses 0, 0|, we replace 
the summation over discrete crystalline neighbor loca- 
tion with a mean-field average over the first shell of 
the pair distribution function of the hardsphere liquids, 
agi{a,n) = % j^,t,,h,g{R.n) TrVVv{R-a)dR. For nu- 
merical work, we take g{R,n) as the Verlet and Weis's 
corrected radial distribution function . 

After some calculations, wc obtain the steady-state 
manybody probability distribution function as a product 
of localized Gaussians of the form exp[— Q!(fi — with 
a, the final localization strength, satisfying two equa- 
tions: 
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of eqn. (5) come from 
respectively. The integral 
^ can be 
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further expressed as explicit but complicated analytical 
formulas with incomplete Gamma functions. Thus using 
eqn. (4) in eqn. (5), we finally derive the order parame- 
ter for localization, a, in a closed form with parameters 
£, D, K, s, s' and n. When k ■ £ ~ Q, wc have ^ = and 
the equilibrium equation a — a(a, n), which returns to 
the self-consistent plionon solution 0, J^, _1£, 2^ ■ A 
nonzero solution a of eqn. (5) is only obtained at suffi- 
ciently high density, i.e., for n > ric- For low density, the 
system cannot support stable localized vibrations and is 
in the fluid phase with 5 = 0. An instability density 
ric separates these two phases. This phase transition is 
first-order like, characterized by a discontinuous jump of 
a. 

We calculated function of two independent pa- 

rameters s and s', nc(s,s') for various k, D, and £. An 
important dimensionless ratio A := k£^ /D measures the 
strength of chemical vs. thermal noise. For an actin 
polymer solution 0], we can relate the effective kick- 
ing rate to the speed of nonequilibrium polymerization. 
Here £ is the monomer size 0.01/.(m. (// = 10^^) The 
treadmill concentration is ctr ~ 0.17/iM. The chem- 
ical reaction rates of the barbed and pointed ends of 



actin are — 11.6/iAf s and k. 
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respectively. The diffusion constants of rod are given 
by the Kirkwood equation This equation relates 
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the diffusion to solvent viscosity rjs and gives the trans- 
lational diffusion constants: Z3|| = kgT ln{L / d) / {2Trr]sL) 
and D± = -D||/2. Using the typical length of the actin fil- 
aments L 20^m and typical width d of about 0.015//to. 
An estimate of the hydrodynamic diffusion along the rod 
gives D|[ ~ Q.lfim^s^^ . Thus for a dilute solution of typ- 
ical actin filaments, A ~ 10~'^-10~^. However, in vivo, 
the actin monomer concentration is kept much higher 
than ctr (with the help of capping proteins which prevent 
actin filaments from growing longer). Also the viscosity 
of the cell medium is higher than rjg of pure water due to 
the presence of other macromolecular components. These 
lower the effective diffusion constant and raise the value 
of K, and therefore they could push A above 1. The limit 
A 3> 1 corresponds to entirely motorized motion. 

The resulting densities nc(s, s') are shown in Fig. 1 and 
2. Figure 1 shows some ID plots that come from vertical 
slices of Uc for several special cases. Figure 2 shows the 
2D top-view of the critical surfaces ric for variety of pa- 
rameters. For the case shown with A S> 1, a{n) (which 
is the only n dependent part of a) drops to zero, while 
we still have a nonzero solution (5 = ^ for a corner of 
(s, s') space. In the opposite corner we stopped searching 
for solutions when Uc approached the maximum packing 
density ^/2. Thus we have three distinct regions. 

As seen from these figures, kicking noise does not al- 
ways destabilize the structures. Instead the localized 
phases have an enlarged stability range when s + s' > 1. 
When s' = 1 — s, the same stability limits arc obtained 
as in the equilibrium thermodynamic theory. Both the 
frozen disordered glass and the ordered fee lattice can be 
stabilized by kicking motors. The motor effects on the fee 
lattice are more pronounced. The fee phase has a larger 
stable region than the glass. 

Besides the Eyink variational method, we also calcu- 
lated a, a, and therefore Uc by another method. From 
the mechanical feedback procedure, we first obtain the 
mapping from a hardsphere environment to an effec- 
tive harmonic potential a = a{a). Here a depends on 
the steady-state probability distribution of its neighbors. 
Conversely a can be viewed as the final effective spring 
constant of a kicking particle in an harmonic potential of 
a. Next we numerically solve a ~ a{a) from single parti- 
cle master equations using a variational method based on 
the square hermitianized operator I'^l with single particle 
trial hmctions. The two sets of operations are iterated to 
obtain a pair of self-consistent results (a, a). The critical 
density predicted by this self-consistent squared hermi- 
tian variational method agrees very well with results from 
the nonhermitian variational method. The difference of 
instability density is less then 0.1% when A < 1. The 
corresponding critical etc are also similar. The two meth- 
ods do give different results when A::S> 1. 

For persistent kicking, the trial functions have to be 
modified. Each localized particle now has a distribution 
of locations of the form ip^{r ' = r — 6; a) ~ exp(— a : 
f'f'). Here b is an off-center shift vector parallel to n. 
We must consider the effects of the variational parameter 
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FIG. 1: The instability density of the motorized fee lattice 
and the glass as functions of the coupling parameters for these 
cases i) s = ii) s' = and iii) s = s' . In these plots there 
is 3D isotropic kicking with Z) = 0.1, k = 10, and I = 0.05. 
Therefore A = = 0.25. The two horizontal lines are the 
corresponding equilibrium [kI — 0) cases. 
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FIG. 2: The instability densities nc{s,s') surface are shown 
for motorized fee lattice cases: (a) A — 0.25, (b) A = 2.5, 
(c), A = 25 and a glass case (d) A — 0.25 for the 3D isotropic 
kicking with k = 10 and £ = 0.05. The corresponding surface 
for the Lindemann parameter ac{s,s') of the case of (c) is 
shown in (e). 
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FIG. 3: The instability density of the persistent motorized 
fee lattice and the glass as parametric functions of h, which 
depends on A. Here I = 0.05. Both are bounded with middle 
line s + s' = 1, with upper and lower bound s = s' = and 1 
respectively. 

6 on a along with the direct changes of a. The additional 
decrease of a arises from the distortion of the structures 
caused by always kicking in the same directions. We 
model the distortion effect of persistent kicking on the 
pair distributions by replacing each initial position with 
a dispersed distribution. For the crystal, this means the 
initial neighboring position of Rj is replaced by an av- 
erage over positions Rj + hn with h is an arbitrary unit 
direction. Likewise, the radial distribution function of 
the glass case is broadened from the initial gh=Q(r) to 
9b{f) = J 9b=o{r'){j^)S{\r' - r\ - b)dr' . In this case 
an additional normalization of the first peak enforces the 
condition g{r) = for r < 1. 

For persistent kicking, s has similar effects on a as 
were found for the isotropic case. The shift b agrees 
quite well (within several percentages for the practical 
range of parameters) with the estimate based on a 
small £ expansion of the master equation. The resulting 
displacement amplitude b is large enough to distort the 
stable structure so that any increased stability that may 
arise from kicking (if any) is now very modest as shown 
in Fig. 3. 



Since the kicking noise enlarges the stability region in 
the isotropic case with s + s' > 1, we wondered whether 
susceptible kicking may sometimes actually decrease the 
effective temperature of this nonequilibrium system. An 
effective temperature can be defined by the fluctuation- 
dissipation relation even for far-from-equilibrium situa- 
tions like the motorized crystal The ratio of the 
thermal temperature to the effective temperature is also 
called FDT violation factor. The effective temperature 
can be computed by comparing the fluctuations of a mo- 
torized particle to its response to an external force. T'^-^^ 
depends on the frequency or time duration, the absolute 
time (in case of an aging system) and even possibly on the 
choice of the observable itself. To compute the needed 
time-dependent quantities, we solve the time dependent 
master equations for nonhcrmitian operators, again us- 
ing a Gaussian ansatz characterized now with dynamic 
first and second moments. 

For illustration, we carry out the analysis of 
the dynamics for the one dimensional symmet- 
ric case in a harmonic potential ax^ with s = 
s'. Thus LpF tjj = Ddl^i) - d^[l3Df''{x)^] 
and Lne i^ix) = fV'(a; - £)e^(^-^/) + |^(a; + 

^ye(x+L-l) _ |^(a;)[es(^,^) + ee{,x.-t)^ ^j^j^ ^ 

~2saxl — saP . The parameters in the Gaussian 
ansatz V^(x,m(t)) = [2^m2(t)]"^/^ exp (- '^'^^gl' ) 

with TO2 := m2 — m\ and nii = {x^) ^ sat- 
isfy the time-dependent dynamics described by a set 
of differential equations with K, = k exp(— sai?^): 
dtmi = -2Dami - /C&^f''"'^)'"^ sinh(2sa^mi) and 
dtm2 ^ 2D - 4Dam2 + O^e^f''"^)'™^ cosh(2sa^mi) - 
20e2('*"'^)'"2 sinh(2sa^TOi) + 2sa^TO2 cosh(2sa^TOi)]. 

To obtain the Green's function G{x, x'; t, 0), mi{t) and 
m2{t) must satisfy the above equations with the initial 
conditions toi(0) = x' and 7712(0) = x' . We denote 
m(t; x') for this pair solution of the differential equa- 
tion. It is easy to see that for k = 0, the exact expres- 
sion of the Green's function of an equilibrium system 
is recovered. Green's function yields the correlations 
and responses. The correlation functions are given by 
C{t) ~ J dx'x'mi{t;x')'4'ix' ,m*). Here * donates the 
steady-state value. While the response to a pulse is 
R{t) = J dx'{l3D/m^)x'mi{t;x')ip{x',in*). Combining 
these yields the effective temperature 



T^ffjt) _ ^ dtC[t) _ ( l + ^^exp[2(sa£)'^m*] ^ ^ ( _^ ^ JC£ J e^(''"^)''^^(*-^') smh{2sa£mi{t; x'))x'ip*{x')dx' 



T*^ R{t) V l + s^exp[2(sa£)2m*] / \^ 2Da J mi[t;x')x'^*{x')dx' ) 



Thus we see that FDT violation is a product of two with the corresponding thermal equilibrium value. The 
ratios. One ratio is the steady-state variance compared other ratio depends on the rate at which the system 
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reaches the steady state. I.e., the larger the variance 
and the faster the dynamics, the hotter the system and 
vice versa. Compared with the cases without kicking, 
susceptible motors yield a smaller variance while on the 
other hand, they relax faster. In the short time limit, 
mi{t;x') — > x' and rh2(t;x') — > 0, and the second ratio 
becomes 1 + 5Mlg2(saf)^ni2 _ Pqj. ^Qjjg times, mi 0, 
sinh(sa^TOi) — > sa£mi, and 7Ji2 — > 7713. We find at this 

jk-/,2 (sat)'' » 

limit the value of the ratio as 1 + s^^e 2 exact 
same as the short time limit. Therefore T'^ffit = 0) = 
T^f-fit = 00) = 1 + 1 «!e2(sQ£)=m*^ Yet for intermediate 
times, the ratio is not constant and differs from either 
limiting value. Generally, T^ff > T*^, i.e., the system 
is "hotter" although chemical noise apparently enlarges 
the stability range of the localized phase. 

We have studied the stability and dynamics of localized 
nonequilibrium structures of motorized particles. The 



nonequilibrium noise from kicking motors sometimes in- 
creases the effective spring constant and enlarges the me- 
chanical stability range of both crystal and frozen glass 
structures. We see that for systems like the cytoskeleton, 
nonequilibrium noise may speed up the dynamics without 
sacrificing structural stability. This model can be further 
developed to include anisotropy of the particles or under 
other types of nonequilibrium noise or driven forces. We 
have also studied the dynamics and effective tempera- 
tures j23| of this system. Besides taking this solid-state 
viewpoint, one can also study the transition from the liq- 
uid side by mode coupling theory, a problem for future 
works. 
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especially in the early stage of the work. 
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